---
title: "Perceptions of the Social Status Hierarchy and Its Cultural and Economic Sources"
subtitle: "Replication File"
author: "Magdalena Breyer"
date: "2024-06-05"
output: 
    html_document:
        toc: true
        code_folding: hide
---

```{r setup, include=TRUE}
knitr::opts_chunk$set(echo = TRUE, warning = F, message = F)

# packages
if (!require("pacman")) install.packages("pacman") 

pacman::p_load(tidyverse, # v2.0.0
               cregg, # v0.4.0
               scales, # v1.3.0
               viridis, # v0.6.5
               RColorBrewer, # v1.1.3
               gridExtra, # v2.3
               kableExtra, # v1.4.0
               modelsummary # v2.1.0
)
theme_set(theme_bw())


```


# Data 

```{r}

# choice task

df_long <- readRDS("../data/df_choice_rep.rds")

# rating task

df_rate_long <- readRDS("../data/df_rating_rep.rds") 

# wide format (for descriptive statistics, NA not excluded)

df_wide <- readRDS("../data/df_wide_rep.rds") 
```


# Main text results

```{r}
# Formula choice and rating outcome

fc <- selected ~  Occupation * Education + Income + Gender + Sexuality +
  Migration.background + Place.of.residence + Age

```


## Figure 1

```{r, fig.align='center', fig.width=6, fig.height=7, fig.cap="Marginal means for the sources of social status. Note: Results of the choice-based outcome of the conjoint. Standard errors were clustered to level of respondents and bars display 95% confidence intervals."}
statusmm <- cj(df_long, estimate = "mm", fc, id = ~ respondent)

plot(statusmm, vline = 0.5) + 
  scale_x_continuous(limits = c(0.25, 0.85), breaks = c(0.3, 0.5, 0.7, 0.9)) +
  theme_bw() +
  theme(legend.position = "none") +
  scale_y_discrete(labels = rev(c("(Occupation)", "Lawyer", "Engineer", "Primary school teacher",
                                  "Electrician", "Gardener", "Waiter",
                                  "(Education)", "University degree", "Vocational degree", "High school degree",
                                  "(Income)", "High income", "Medium income", "Low income",
                                  "(Gender)", "Female", "Male",
                                  "(Sexuality)", "Homosexual", "Heterosexual",
                                  "(Migration background)", "Nigerian", "Kosovar", "Turkish", "Italian", "None",
                                  "(Place of residence)", "Countryside", "Suburbs", "Big city",
                                  "(Age)", "62 years", "46 years", "29 years")))

```


## Figure 2

```{r, fig.align='center', fig.width=6, fig.height=5, fig.cap="Importance of attributes by status dimension. Note: Importance in percent denotes the share an attribute contributes to the range of AMCEs of all attributes"}

# AMCEs

statusamces <- cj(df_long, fc, estimate = "amce",
                  id = ~ respondent) %>% 
  mutate(
    dimension = if_else(feature %in% c("Occupation", "Education", "Income"), "Economic",
                        if_else(feature %in% c("Gender", "Sexuality", "Migration.background", "Place.of.residence"), "Cultural",
                                "Neither")),
    dimension = relevel(factor(dimension), "Economic", "Cultural", "Neither")
  )


# take into account bounds to standardize AMCEs

statusamces_s <- statusamces %>%
  group_by(feature) %>%
  add_tally() %>% # number of levels per feature
  mutate(
    low_bound = - (1 - (1 / n)),
    up_bound = (1 - (1 / n)),
    estimate_old = estimate,
    estimate = rescale(estimate, to = c(-1, 1), from = range(low_bound, up_bound)),
    lower = rescale(lower, to = c(-1, 1), from = range(low_bound, up_bound)), #CIs
    upper = rescale(upper, to = c(-1, 1), from = range(low_bound, up_bound))
    ) 

# Range of AMCEs in one attribute as percentages of total range 

feat_imp_s <- statusamces_s %>%
  group_by(feature) %>%
  summarise(maxe = max(estimate),
            mine = min(estimate),
            rng = maxe-mine) %>%
  ungroup %>%
  mutate(
    sumr = sum(rng),
    importance = (rng / sumr) * 100,
    sumi = sum(importance)
  ) %>%
  mutate(
    dimension = if_else(feature %in% c("Occupation", "Education", "Income"), "economic",
                        if_else(feature %in% c("Gender", "Sexuality", "Migration.background", "Place.of.residence"), "cultural",
                                "neither"))
  )

ggplot(feat_imp_s, aes(fct_reorder(feature, importance, .desc = F), importance, fill = dimension)) + 
  geom_col() + coord_flip() +
  scale_fill_grey() +
  scale_y_continuous(limits = c(0,30)) +
  theme(legend.position = "bottom") +
  theme(legend.title=element_blank()) +
  labs(x = NULL, y = "Importance in %") +
  scale_x_discrete(labels = rev(c("Occupation", "Migration background", "Income", "Education", 
                              "Sexuality", "Gender", "Age", "Place of residence")))

```



## Figure 3

```{r, fig.width=9, fig.height=7, fig.cap="Subgroup heterogeneity by respondent migration background. Note: Results of the choice-based outcome of the conjoint: Left plot shows marginal means by respondent migration background, right plot shows difference between the two groups. Standard errors were clustered to level of respondents and bars display 95% confidence intervals."}

mm_by_mig <- cj(df_long, fc, estimate = "mm", id = ~ respondent,
                by = ~ migrationb)

mm_by_mig_diff <- cj(df_long, fc, estimate = "mm_diff", id = ~ respondent,
                     by = ~ migrationb)

smig <- plot(mm_by_mig, group = "migrationb", vline = 0.5, size = 1) +
  scale_color_brewer(palette = "Paired", na.translate = F) +
  theme_bw() +
  geom_hline(yintercept = c(4, 8, 14, 17, 20, 24, 28), size = 0.25) +
  scale_y_discrete(labels = rev(c("(Occupation)", "Lawyer", "Engineer", "Primary school teacher",
                                  "Electrician", "Gardener", "Waiter",
                                  "(Education)", "University degree", "Vocational degree", "High school degree",
                                  "(Income)", "High income", "Medium income", "Low income",
                                  "(Gender)", "Female", "Male",
                                  "(Sexuality)", "Homosexual", "Heterosexual",
                                  "(Migration background)", "Nigerian", "Kosovar", "Turkish", "Italian", "None",
                                  "(Place of residence)", "Countryside", "Suburbs", "Big city",
                                  "(Age)", "62 years", "46 years", "29 years"))) +
  theme(legend.position = "bottom", legend.title=element_blank(),
        text = element_text(size=10)) +
  guides(color=guide_legend(ncol=2,nrow=1,byrow=T))


smig_diff <- plot(mm_by_mig_diff, group = "migrationb", size = 1) +
  geom_hline(yintercept = c(4, 8, 14, 17, 20, 24, 28), size = 0.25) +
  scale_y_discrete(labels = NULL) +
  scale_colour_grey(na.translate = F) + 
  theme_bw() +
  theme(legend.position = "bottom", legend.title=element_blank()) +
  theme(text = element_text(size=10))


pmig <- arrangeGrob(smig, smig_diff, ncol = 2,
                     widths = unit(c(10,7), c("cm")))
grid.arrange(pmig)
```

## Figure 4

```{r, fig.width=9, fig.height=7, fig.cap="Subgroup heterogeneity by respondent urban or rural residence. Note: Results of the choice-based outcome of the conjoint: Left plot shows marginal means by respondent place of residence, right plot shows difference between the two groups. Standard errors were clustered to level of respondents and bars display 95% confidence intervals."}

mm_by_ur <- cj(df_long, fc, estimate = "mm", id = ~ respondent,
               by = ~ urban_2cat)

mm_by_ur_diff <- cj(df_long, fc, estimate = "mm_diff", id = ~ respondent,
                    by = ~ urban_2cat)

sur <- plot(mm_by_ur, group = "urban_2cat", vline = 0.5, size = 1) +
  scale_color_brewer(palette = "Paired", na.translate = F) +
  theme_bw() +
  geom_hline(yintercept = c(4, 8, 14, 17, 20, 24, 28), size = 0.25) +
  scale_y_discrete(labels = rev(c("(Occupation)", "Lawyer", "Engineer", "Primary school teacher",
                                  "Electrician", "Gardener", "Waiter",
                                  "(Education)", "University degree", "Vocational degree", "High school degree",
                                  "(Income)", "High income", "Medium income", "Low income",
                                  "(Gender)", "Female", "Male",
                                  "(Sexuality)", "Homosexual", "Heterosexual",
                                  "(Migration background)", "Nigerian", "Kosovar", "Turkish", "Italian", "None",
                                  "(Place of residence)", "Countryside", "Suburbs", "Big city",
                                  "(Age)", "62 years", "46 years", "29 years"))) +
  theme(legend.position = "bottom", legend.title=element_blank(),
        text = element_text(size=10))

sur_diff <- plot(mm_by_ur_diff, group = "urban_2cat", size = 1) +
  geom_hline(yintercept = c(4, 8, 14, 17, 20, 24, 28), size = 0.25) +
  scale_y_discrete(labels = NULL) +
  scale_colour_grey(na.translate = F) + 
  theme_bw() +
  theme(legend.position = "bottom", legend.title=element_blank()) +
  theme(text = element_text(size=10))

psur <- arrangeGrob(sur, sur_diff, ncol = 2,
                     widths = unit(c(10,7), c("cm")))
grid.arrange(psur)

```


# Online Appendix

## Table A1

```{r}
# categorical
datasummary(data = df_wide,
            `Gender` + `Income group (household equivalized)` + 
              `Education group` + `Age group` +
              `Sexual orientation` + `Migration background` + `Urban/rural residence` +
              `Social status group` + `Past social status group` + `Future social status group` +
              `Social dominance orientation group` + `Left/right self-placement groups` + 1 ~  N + Percent(),
            title = "Sample descriptives, categorical variables",
            output = "kableExtra"
            ) %>% kable_classic()
```

## Table A2

```{r}
# continuous
datasummary_skim(data = df_wide, type = "numeric",
                 title = "Sample descriptives, categorical variables",
                 output = "kableExtra"
                 ) %>% kable_classic()
```

## Figure D1

```{r, fig.align='center', fig.width=6, fig.height=7, fig.cap="Average marginal component effects for the sources of social status, choice-based outcome."}

plot(statusamces, size = 1) +
  theme_bw() +
  theme(legend.position = "none") +
  scale_y_discrete(labels = rev(c("(Occupation)", "Lawyer", "Engineer", "Primary school teacher",
                              "Electrician", "Gardener", "Waiter",
                              "(Education)", "University degree", "Vocational degree", "High school degree",
                              "(Income)", "High income", "Medium income", "Low income",
                              "(Gender)", "Female", "Male",
                              "(Sexuality)", "Homosexual", "Heterosexual",
                              "(Migration background)", "Nigerian", "Kosovar", "Turkish", "Italian", "None",
                              "(Place of residence)", "Countryside", "Suburbs", "Big city",
                              "(Age)", "62 years", "46 years", "29 years")))
```

## Table D1

```{r}

statusamces %>% 
  mutate(feature = fct_recode(feature, 
                              "Migr.backg" = "Migration.background",
                              "Place" = "Place.of.residence")) %>% 
  select("Attribute" = feature, "Level" = level, 
         "AMCE" = estimate, "Std. Error" = std.error, p) %>% 
  mutate(across(where(is.numeric), ~as.character(signif(., 3)))) %>% 
  kbl(digits = 3, caption = "Choice-based outcome, AMCEs") %>% 
  kable_classic(full_width = F)
```

## Table D2

```{r}
statusmm %>% 
  mutate(feature = fct_recode(feature,
                              "Migr.backg" = "Migration.background",
                              "Place" = "Place.of.residence")) %>%
  select("Attribute" = feature, "Level" = level, 
         "Marginal mean" = estimate, "Std. error" = std.error, p) %>% 
  mutate(across(where(is.numeric), ~as.character(signif(., 3)))) %>% 
  kbl(digits = 3, caption = "Choice-based outcome, marginal means") %>% 
  kable_classic(full_width = F)

```


## Table E1

```{r}
statusamces_s %>% 
  mutate(
    Bounds = paste0(round(low_bound, 2), ";", round(up_bound, 2))
  ) %>% 
  select("Attribute" = feature, "Level" = level, "# Levels" = n, Bounds,
         "AMCE" = estimate_old, "rescaled AMCE" = estimate) %>% 
  mutate(across(where(is.numeric), ~as.character(signif(., 3)))) %>% 
  kbl(digits = 3, caption = "Choice-based outcome, rescaled AMCEs") %>% 
  kable_classic(full_width = F)
```


## Table E2

```{r}
feat_imp_s %>% 
  mutate(feature = fct_recode(feature,
                              "Migr.backg" = "Migration.background",
                              "Place" = "Place.of.residence")) %>%
  arrange(-importance) %>% 
  select("Attribute" = feature, "Dimension" = dimension, "AMCE range" = rng,  
         "Sum of ranges" = sumr, "Importance in %" = importance) %>% 
  kbl(digits = 3, caption = "Choice-based outcome, importance of attributes") %>% 
  kable_classic(full_width = F)
```

## Figure F1

```{r, fig.align='center', fig.width=6, fig.height=7, fig.cap="Rating outcome: Marginal Means"}


rstatusmm <- cj(df_rate_long, estimate = "mm", fc, id = ~ respondent)

plot(rstatusmm, vline = 6.48) +
  scale_x_continuous(breaks = c(5.5, 6.5, 7.5, 8.5)) +
  theme_bw() +
  theme(legend.position = "none") +
  scale_y_discrete(labels = rev(c("(Occupation)", "Lawyer", "Engineer", "Primary school teacher",
                                  "Electrician", "Gardener", "Waiter",
                                  "(Education)", "University degree", "Vocational degree", "High school degree",
                                  "(Income)", "High income", "Medium income", "Low income",
                                  "(Gender)", "Female", "Male",
                                  "(Sexuality)", "Homosexual", "Heterosexual",
                                  "(Migration background)", "Nigerian", "Kosovar", "Turkish", "Italian", "None",
                                  "(Place of residence)", "Countryside", "Suburbs", "Big city",
                                  "(Age)", "62 years", "46 years", "29 years")))
```

## Figure F3

```{r, fig.align='center', fig.width=6, fig.height=7, fig.cap="Rating outcome: AMCEs"}
rstatusamces <- cj(df_rate_long, fc, estimate = "amce",
                  id = ~ respondent)


plot(rstatusamces, size = 1) +
  theme_bw() +
  theme(legend.position = "none") +
  scale_y_discrete(labels = rev(c("(Occupation)", "Lawyer", "Engineer", "Primary school teacher",
                                  "Electrician", "Gardener", "Waiter",
                                  "(Education)", "University degree", "Vocational degree", "High school degree",
                                  "(Income)", "High income", "Medium income", "Low income",
                                  "(Gender)", "Female", "Male",
                                  "(Sexuality)", "Homosexual", "Heterosexual",
                                  "(Migration background)", "Nigerian", "Kosovar", "Turkish", "Italian", "None",
                                  "(Place of residence)", "Countryside", "Suburbs", "Big city",
                                  "(Age)", "62 years", "46 years", "29 years")))
```

## Table F3

```{r}
rstatusamces %>% 
  mutate(feature = fct_recode(feature, 
                              "Migr.backg" = "Migration.background",
                              "Place" = "Place.of.residence")) %>% 
  select("Attribute" = feature, "Level" = level, 
         "AMCE" = estimate, "Std. Error" = std.error, p) %>% 
  mutate(across(where(is.numeric), ~as.character(signif(., 3)))) %>% 
  kbl(digits = 3, caption = "Rating outcome, AMCEs") %>% 
  kable_classic(full_width = F)
```

```{r}
rstatusmm %>% 
  mutate(feature = fct_recode(feature, 
                              "Migr.backg" = "Migration.background",
                              "Place" = "Place.of.residence")) %>% 
  select("Attribute" = feature, "Level" = level, 
         "Marginal mean" = estimate, "Std. error" = std.error, p) %>% 
  mutate(across(where(is.numeric), ~as.character(signif(., 3)))) %>% 
  kbl(digits = 3, caption = "Rating outcome, marginal means") %>% 
  kable_classic(full_width = F)
```

## Figure G1

```{r, fig.align='center', fig.width=6, fig.height=7, fig.cap="Interaction effects of profile migration background with all other profile attributes, marginal means, choice outcome"}

statusmm_pi <- cj(df_long, estimate = "mm", id = ~ respondent,
                  selected ~ Occupation * Education + Income + Gender + Sexuality + 
                    Place.of.residence + Age,
                  by = ~Migration.background)


plot(statusmm_pi, group = "Migration.background", vline = 0.5) +
  scale_colour_manual(name = "", values = brewer.pal(9,"Blues")[c(3,5,6,8,9)],
                      na.translate = F) +
  theme_bw() + theme(legend.position = "bottom") +
  theme(legend.text=element_text(size=8))+
  guides(color=guide_legend(ncol=3,nrow=2,byrow=T)) +
  scale_y_discrete(labels = rev(c("(Occupation)", "Lawyer", "Engineer", "Primary school teacher",
                                  "Electrician", "Gardener", "Waiter",
                                  "(Education)", "University degree", "Vocational degree", "High school degree",
                                  "(Income)", "High income", "Medium income", "Low income",
                                  "(Gender)", "Female", "Male",
                                  "(Sexuality)", "Homosexual", "Heterosexual",
                                  "(Place of residence)", "Countryside", "Suburbs", "Big city",
                                  "(Age)", "62 years", "46 years", "29 years")))
```

## Figure G2

```{r, fig.align='center', fig.width=6, fig.height=7, fig.cap="Interaction effects of profile gender with all other profile attributes, marginal means, choice outcome"}

statusmm_pig <- cj(df_long, estimate = "mm", id = ~ respondent,
                   selected ~ Occupation * Education + Income + Sexuality + 
                     Migration.background + Place.of.residence + Age, by = ~Gender)


plot(statusmm_pig, group = "Gender", vline = 0.5, size = 1) +
  scale_color_brewer(palette = "Paired", na.translate = F) +
  guides(shape = "none") + theme_bw() +
  theme(legend.text=element_text(size=8))+
  theme(legend.position = "bottom", legend.title=element_blank())+
  scale_y_discrete(labels = rev(c("(Occupation)", "Lawyer", "Engineer", "Primary school teacher",
                                  "Electrician", "Gardener", "Waiter",
                                  "(Education)", "University degree", "Vocational degree", "High school degree",
                                  "(Income)", "High income", "Medium income", "Low income",
                                  "(Sexuality)", "Homosexual", "Heterosexual",
                                  "(Migration background)", "Nigerian", "Kosovar", "Turkish", "Italian", "None",
                                  "(Place of residence)", "Countryside", "Suburbs", "Big city",
                                  "(Age)", "62 years", "46 years", "29 years")))
```
## Figure G3

```{r, fig.align='center', fig.width=6, fig.height=7, fig.cap="Interaction effects of profile sexual orientation with all other profile attributes, marginal means, choice outcome"}

statusmm_pis <- cj(df_long, estimate = "mm", id = ~ respondent,
                   selected ~ Occupation * Education + Income + Gender + 
                     Migration.background + Place.of.residence + Age, by = ~Sexuality)


plot(statusmm_pis, group = "Sexuality", vline = 0.5, size = 1) +
  scale_color_brewer(palette = "Paired", na.translate = F) +
  guides(shape = "none") + theme_bw() +
  theme(legend.text=element_text(size=8))+
  theme(legend.position = "bottom", legend.title=element_blank())+
  scale_y_discrete(labels = rev(c("(Occupation)", "Lawyer", "Engineer", "Primary school teacher",
                                  "Electrician", "Gardener", "Waiter",
                                  "(Education)", "University degree", "Vocational degree", "High school degree",
                                  "(Income)", "High income", "Medium income", "Low income",
                                  "(Gender)", "Female", "Male",
                                  "(Migration background)", "Nigerian", "Kosovar", "Turkish", "Italian", "None",
                                  "(Place of residence)", "Countryside", "Suburbs", "Big city",
                                  "(Age)", "62 years", "46 years", "29 years")))
```

## Figure G4

```{r, fig.align='center', fig.width=6, fig.height=7, fig.cap="Interaction effects of profile place of residence with all other profile attributes, marginal means, choice outcome"}

# place
statusmm_pip <- cj(df_long, estimate = "mm", id = ~ respondent,
                   selected ~ Occupation * Education + Income + Gender + Sexuality +
                     Migration.background +  Age, by = ~Place.of.residence)


plot(statusmm_pip, group = "Place.of.residence", vline = 0.5, size = 1) +
  scale_color_manual(values = c("dodgerblue4", "lightskyblue3", "lightcyan2"), na.translate = F) +
  guides(shape = "none") + theme_bw() +
  theme(legend.text=element_text(size=8))+
  theme(legend.position = "bottom", legend.title=element_blank())+
  scale_y_discrete(labels = rev(c("(Occupation)", "Lawyer", "Engineer", "Primary school teacher",
                                  "Electrician", "Gardener", "Waiter",
                                  "(Education)", "University degree", "Vocational degree", "High school degree",
                                  "(Income)", "High income", "Medium income", "Low income",
                                  "(Gender)", "Female", "Male",
                                  "(Sexuality)", "Homosexual", "Heterosexual",
                                  "(Migration background)", "Nigerian", "Kosovar", "Turkish", "Italian", "None",
                                  "(Age)", "62 years", "46 years", "29 years")))

```

## Table G1

```{r}

statusmm_pi %>% 
  mutate(feature = fct_recode(feature, 
                              "Place" = "Place.of.residence")) %>% 
  select("By PMB" = BY, "Attribute" = feature, "Level" = level, 
         "MM" = estimate, "Std. error" = std.error, p) %>% 
  mutate(across(where(is.numeric), ~as.character(signif(., 3)))) %>% 
  kbl(digits = 3, caption = "Interaction effects of profile migration background (PMB) with all other attributes, marginal means (MM)") %>% 
  kable_classic(full_width = F)


```

## Table H1

```{r}

# ANOVA per respondent subgroup
anovasss <- df_long %>%
  filter(!is.na(sss_group2)) %>%
  cj_anova(fc, by = ~sss_group2) %>%
  mutate(subgroup = "Subjective social status")

anovag <- df_long %>%
  filter(!is.na(gender)) %>%
  cj_anova(fc, by = ~gender) %>%
  mutate(subgroup = "Gender")

anovas <- df_long %>%
  filter(!is.na(sexual_orientation)) %>%
  cj_anova(fc, by = ~sexual_orientation) %>%
  mutate(subgroup = "Sexual orientation")

anovam <- df_long %>%
  filter(!is.na(migrationb)) %>%
  cj_anova(fc, by = ~migrationb) %>%
  mutate(subgroup = "Migration background")

anovai <- df_long %>%
  filter(!is.na(income_group_eq)) %>%
  cj_anova(fc, by = ~income_group_eq) %>%
  mutate(subgroup = "Income")

anovae <- df_long %>%
  filter(!is.na(edu)) %>%
  cj_anova(fc, by = ~edu) %>%
  mutate(subgroup = "Education")

anovaa <- df_long %>%
  filter(!is.na(age_group)) %>%
  cj_anova(fc, by = ~age_group) %>%
  mutate(subgroup = "Age")

anovaur <- df_long %>%
  filter(!is.na(urban_2cat)) %>%
  cj_anova(fc, by = ~urban_2cat) %>%
  mutate(subgroup = "Urban/rural residence")

anovasdo <- df_long %>%
  filter(!is.na(sdo_group)) %>%
  cj_anova(fc, by = ~sdo_group) %>%
  mutate(subgroup = "Social dominance drientation")

anovalr <- df_long %>%
  filter(!is.na(lr_group)) %>%
  cj_anova(fc, by = ~lr_group) %>%
  mutate(subgroup = "Left/right self-positioning")


anova_all <- bind_rows(anovasss, anovasdo, anovalr, anovai, anovas, anovam, anovag, anovae, 
                       anovaa, anovaur)  %>%
  mutate_if(is.numeric, format, digits=2)

anova_all %>% 
  arrange(subgroup) %>%
  select(-subgroup) %>%
    kbl(caption = "Analysis of deviance: Test for subgroup heterogeneity based on respondent
characteristics (economic and cultural sources)"
    ) %>% 
  kable_classic(full_width = F) %>%
  pack_rows(index = table(anova_all$subgroup))
```

## Figure H1


```{r, fig.width=9, fig.height=7, fig.cap="Subgroup heterogeneity by respondent subjective social status. Note: Results of the choice-based outcome of the conjoint: Left plot shows marginal means by respondent SSS, right plot shows difference between the two groups. SSS was surveyed on a scale from 1 to 10, those with status below the median (7) are low status, while those with values of 7 or above are high status."}

mm_by_sss <- cj(df_long, fc, estimate = "mm", id = ~ respondent,
                by = ~ sss_group2)

mm_by_sss_diff <- cj(df_long, fc, estimate = "mm_diff", id = ~ respondent,
                     by = ~ sss_group2)

ssss <- plot(mm_by_sss, group = "sss_group2", vline = 0.5, size = 1) +
  scale_color_brewer(palette = "Paired", na.translate = F) +
  theme_bw() +
  geom_hline(yintercept = c(4, 8, 14, 17, 20, 24, 28), size = 0.25) +
  scale_y_discrete(labels = rev(c("(Occupation)", "Lawyer", "Engineer", "Primary school teacher",
                                  "Electrician", "Gardener", "Waiter",
                                  "(Education)", "University degree", "Vocational degree", "High school degree",
                                  "(Income)", "High income", "Medium income", "Low income",
                                  "(Gender)", "Female", "Male",
                                  "(Sexuality)", "Homosexual", "Heterosexual",
                                  "(Migration background)", "Nigerian", "Kosovar", "Turkish", "Italian", "None",
                                  "(Place of residence)", "Countryside", "Suburbs", "Big city",
                                  "(Age)", "62 years", "46 years", "29 years"))) +
  theme(legend.position = "bottom", legend.title=element_blank()) +
  theme(text = element_text(size=10))

ssss_diff <- plot(mm_by_sss_diff, group = "sss_group2", size = 1) +
  geom_hline(yintercept = c(4, 8, 14, 17, 20, 24, 28), size = 0.25) +
  scale_y_discrete(labels = NULL) +
  scale_colour_grey(na.translate = F) + 
  theme_bw() +
  theme(legend.position = "bottom", legend.title=element_blank()) +
  theme(text = element_text(size=10))

hsss <- arrangeGrob(ssss, ssss_diff, ncol = 2,
                    widths = unit(c(10,7), c("cm")))
grid.arrange(hsss)
```

## Figure H2



```{r, fig.width=9, fig.height=7, fig.cap="Subgroup heterogeneity by respondent left/right self-placement. Note: Results of the choice-based outcome of the conjoint: Left plot shows marginal means by respondent left/right ideological self-placement, right plot shows difference between the two groups. Left/right self-placement was surveyed on a scale from 0 (left) to 10 (right). The left group comprises all respondents with a value below the median (5), the right group ranges from 5 to 10"}

mm_by_lr <- cj(df_long, fc, estimate = "mm", id = ~ respondent,
                by = ~ lr_group)

mm_by_lr_diff <- cj(df_long, fc, estimate = "mm_diff", id = ~ respondent,
                     by = ~ lr_group)


slr <- plot(mm_by_lr, group = "lr_group", vline = 0.5, size = 1) +
  scale_color_brewer(palette = "Paired", na.translate = F, labels = c("left", "right")) +
  theme_bw() +
  geom_hline(yintercept = c(4, 8, 14, 17, 20, 24, 28), linewidth = 0.25) +
  scale_y_discrete(labels = rev(c("(Occupation)", "Lawyer", "Engineer", "Primary school teacher",
                                  "Electrician", "Gardener", "Waiter",
                                  "(Education)", "University degree", "Vocational degree", "High school degree",
                                  "(Income)", "High income", "Medium income", "Low income",
                                  "(Gender)", "Female", "Male",
                                  "(Sexuality)", "Homosexual", "Heterosexual",
                                  "(Migration background)", "Nigerian", "Kosovar", "Turkish", "Italian", "None",
                                  "(Place of residence)", "Countryside", "Suburbs", "Big city",
                                  "(Age)", "62 years", "46 years", "29 years"))) +
  theme(legend.position = "bottom", legend.title=element_blank()) +
  theme(text = element_text(size=10))

slr_diff <- plot(mm_by_lr_diff, group = "lr_group", size = 1) +
  geom_hline(yintercept = c(4, 8, 14, 17, 20, 24, 28), linewidth = 0.25) +
  scale_y_discrete(labels = NULL) +
  scale_colour_grey(na.translate = F, labels = c("right")) + 
  theme_bw() +
  theme(legend.position = "bottom", legend.title=element_blank()) +
  theme(text = element_text(size=10))

hlr <- arrangeGrob(slr, slr_diff, ncol = 2,
                    widths = unit(c(10,7), c("cm")))
grid.arrange(hlr)
```

## Figure H3

```{r, fig.width=9, fig.height=7, fig.cap="Subgroup heterogeneity by respondent gender. Note: Results of the choice-based outcome of the conjoint: Left plot shows marginal means by respondent gender, right plot shows difference between the two."}

mm_by_gender <- cj(df_long, fc, estimate = "mm", id = ~ respondent,
                   by = ~ gender)

mm_by_gender_diff <- cj(df_long, fc, estimate = "mm_diff", id = ~ respondent,
                        by = ~ gender)

sg <- plot(mm_by_gender, group = "gender", vline = 0.5, size = 1) +
  geom_hline(yintercept = c(4, 8, 14, 17, 20, 24, 28), size = 0.25) +
  scale_color_brewer(palette = "Paired", na.translate = F) +
  guides(shape = "none") + theme_bw() +
  theme(legend.position = "bottom", legend.title=element_blank()) +
  theme(text = element_text(size=10))

sg_diff <- plot(mm_by_gender_diff, group = "gender", size = 1) +
  scale_y_discrete(labels = NULL) +
  geom_hline(yintercept = c(4, 8, 14, 17, 20, 24, 28), size = 0.25) +
  scale_colour_grey(na.translate = F) + 
  theme_bw() +
  theme(legend.position = "bottom", legend.title=element_blank()) +
  theme(text = element_text(size=10))

hg <- arrangeGrob(sg, sg_diff, ncol = 2,
                  widths = unit(c(10,7), c("cm")))
grid.arrange(hg)

```

## Figure H4

```{r, fig.width=9, fig.height=7, fig.cap="Subgroup heterogeneity by respondent sexual orientation. Note: Results of the choice-based outcome of the conjoint: Left plot shows marginal means by respondent sexuality, right plot shows difference between the two."}

mm_by_sexuality <- cj(df_long, fc, estimate = "mm", id = ~ respondent,
                      by = ~ sexual_orientation)

mm_by_sexuality_diff <- cj(df_long, fc, estimate = "mm_diff", id = ~ respondent,
                           by = ~ sexual_orientation)

ssex <- plot(mm_by_sexuality, group = "sexual_orientation", vline = 0.5, size = 1) +
  scale_color_brewer(palette = "Paired", na.translate = F) +
  theme_bw() +
  geom_hline(yintercept = c(4, 8, 14, 17, 20, 24, 28), size = 0.25) +
  scale_y_discrete(labels = rev(c("(Occupation)", "Lawyer", "Engineer", "Primary school teacher",
                                  "Electrician", "Gardener", "Waiter",
                                  "(Education)", "University degree", "Vocational degree", "High school degree",
                                  "(Income)", "High income", "Medium income", "Low income",
                                  "(Gender)", "Female", "Male",
                                  "(Sexuality)", "Homosexual", "Heterosexual",
                                  "(Migration background)", "Nigerian", "Kosovar", "Turkish", "Italian", "None",
                                  "(Place of residence)", "Countryside", "Suburbs", "Big city",
                                  "(Age)", "62 years", "46 years", "29 years"))) +
  theme(legend.position = "bottom", legend.title=element_blank(),
        text = element_text(size=10))

ssex_diff <- plot(mm_by_sexuality_diff, group = "sexual_orientation", size = 1) +
  geom_hline(yintercept = c(4, 8, 14, 17, 20, 24, 28), size = 0.25) +
  scale_y_discrete(labels = NULL) +
  scale_colour_grey(na.translate = F) + 
  theme_bw() +
  theme(legend.position = "bottom", legend.title=element_blank()) +
  theme(text = element_text(size=10))

pssex <- arrangeGrob(ssex, ssex_diff, ncol = 2,
                     widths = unit(c(10,7), c("cm")))
grid.arrange(pssex)

```

## Figure H5

```{r, fig.width=9, fig.height=7, fig.cap="Subgroup heterogeneity by respondent education level. Note: Results of the choice-based outcome of the conjoint: Left plot shows marginal means by respondent edcation level, right plot shows difference between medium/high and low education."}

mm_by_edu <- cj(df_long, fc, estimate = "mm", id = ~ respondent,
                by = ~ edu)

mm_by_edu_diff <- cj(df_long, fc, estimate = "mm_diff", id = ~ respondent,
                     by = ~ edu)

se <- plot(mm_by_edu, group = "edu", vline = 0.5, size = 1) +
  geom_hline(yintercept = c(4, 8, 14, 17, 20, 24, 28), size = 0.25) +
  scale_color_manual(values = c("dodgerblue4", "lightskyblue3", "lightcyan2"), na.translate = F) +
  guides(shape = "none") + theme_bw() +
  theme(legend.position = "bottom", legend.title=element_blank()) +
  theme(text = element_text(size=10))

se_diff <- plot(mm_by_edu_diff, group = "edu", size = 1) +
  geom_hline(yintercept = c(4, 8, 14, 17, 20, 24, 28), size = 0.25) +
  scale_colour_grey(na.translate = F) + 
  scale_y_discrete(labels = NULL) +
  theme_bw() +
  theme(legend.position = "bottom", legend.title=element_blank()) +
  theme(text = element_text(size=10))

he <- arrangeGrob(se, se_diff, ncol = 2,
                  widths = unit(c(10,7), c("cm")))
grid.arrange(he)
```

## Figure H6

```{r, fig.width=9, fig.height=7, fig.cap="Subgroup heterogeneity by respondent age group. Note: Results of the choice-based outcome of the conjoint: Left plot shows marginal means by respondent age group, right plot shows difference between middle-aged/older and younger respondents."}

mm_by_age <- cj(df_long, fc, estimate = "mm", id = ~ respondent,
                by = ~ age_group)

mm_by_age_diff <- cj(df_long, fc, estimate = "mm_diff", id = ~ respondent,
                     by = ~ age_group)

sag <- plot(mm_by_age, group = "age_group", vline = 0.5, size = 1) +
  geom_hline(yintercept = c(4, 8, 14, 17, 20, 24, 28), size = 0.25) +
  scale_color_manual(values = c("dodgerblue4", "lightskyblue3", "lightcyan2"), na.translate = F) +
  theme_bw() +
  theme(legend.position = "bottom", legend.title=element_blank()) +
  theme(text = element_text(size=10))

sag_diff <- plot(mm_by_age_diff, group = "age_group", size = 1) +
  geom_hline(yintercept = c(4, 8, 14, 17, 20, 24, 28), size = 0.25) +
  scale_colour_grey(na.translate = F) + 
  scale_y_discrete(labels = NULL) +
  theme_bw() +
  theme(legend.position = "bottom", legend.title=element_blank()) +
  theme(text = element_text(size=10))

hag <- arrangeGrob(sag, sag_diff, ncol = 2,
                   widths = unit(c(10,7), c("cm")))
grid.arrange(hag)
```

## Table H2

```{r}

tmmm <- left_join(mm_by_mig, mm_by_mig_diff,
                  by = c("feature", "level")) %>% 
  mutate(
    BY.x = fct_recode(factor(BY.x), 
                      "no MB" = "no migration background",
                      "with MB" = "with migration background")
  )


tmmm %>% 
  mutate(feature = fct_recode(feature, 
                              "Migr.backg" = "Migration.background",
                              "Place" = "Place.of.residence")) %>% 
  arrange(feature, level) %>% 
  select("By RMB" = BY.x, "Attribute" = feature, "Level" = level, 
         "MM" = estimate.x,  "MM diff." = estimate.y, "Std. error diff." = std.error.y) %>% 
  mutate(across(where(is.numeric), ~as.character(signif(., 3)))) %>% 
  kbl(digits = 3, caption = "Subgroup marginal means (MM) and difference in MM by respondent migration background (RMB)") %>% 
  kable_classic(full_width = F)

```

## Table H3

```{r}

tmmp <- left_join(mm_by_ur, mm_by_ur_diff,
                  by = c("feature", "level"))

tmmp %>% 
  mutate(feature = fct_recode(feature, 
                              "Migr.backg" = "Migration.background",
                              "Place" = "Place.of.residence")) %>% 
  arrange(feature, level) %>% 
  select("By RP" = BY.x, "Attribute" = feature, "Level" = level, 
         "MM" = estimate.x,  "MM diff." = estimate.y, "Std. error diff." = std.error.y) %>% 
  mutate(across(where(is.numeric), ~as.character(signif(., 3)))) %>% 
  kbl(digits = 3, caption = "Subgroup marginal means (MM) and difference in MM by respondent place of residence (RP)") %>% 
  kable_classic(full_width = F)
```

